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RESUMEN 

Se presenta el marco teorico de las diferencias finitas, orientadas al 
calculo de la matriz de rigidez de elementos de section variable; concretamente de 
vigas con dos cartelas, columnas con variation lineal de las dimensiones de la 
section transversal y de disipadores de energia TADAS. 

Los resultados obtenidos se comparan con los que se hallan al aplicar el 
modelo de las dovelas, por un lado, y por otro con los que reportan las tablas de 
Guldan y de Tena Colunga para las vigas acarteladas. 

Se presenta ademas los programas que han sido elaborados en MATLAB, 
para el efecto, de esta manera el lector complementary sus conocimientos. 

Palabras Clave: Diferencias finitas, Operadores en diferencias; Metodo de 
las dovelas. 

ABSTRACT 



The theoretical framework of finite differences , aimed at calculating the 
stiffness matrix elements of variable section presents ; namely two brackets beams 
, columns with linear variation of cross-sectional dimensions and energy dissipators 
TADAS . 

The results obtained are compared with those found in applying the 
model of the segments, on the one hand and with the other tables reporting and 
Tena Guldan Colunga for gusset beams. 
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Programs that have been developed in MATLAB, to the effect, so the 
reader will complement their knowledge is also presented. 

Key words: Finite difference, Difference operators, Method of the segments. 



INTRODUCCION 



1.1. Matriz de Rigidez 

Por concepto se sabe que un termino cualquiera de la matriz de rigidez 
ktj, es el valor de la carga generalizada Q t correspondiente a la deformada 
elemental qj = 1 y las demas nulas. Si se quiere calcular los elementos de la 
primera columna de la matriz de rigidez, se debe calcular el vector de cargas 
generalizadas que corresponde al estado de desplazamiento elemental 
1i = ly It = 0,i 1. De igual manera se realizan las demas columnas. Aguiar 
(2014) 

En esta investigacion se parte de una viga doblemente acartelada con el 
siguiente sistema P-p de dos grados de libertad, un giro en cada extremo como 
indica la figura 1 .Ademas vale aclarar que cada procedimiento a realizar en la viga 
acartelada sera similar para la columna de seccion variable y el aislador sismico 
tipo TADAS como indica en la figura 2. 




Figura 1 Sistema P-p 




L 



Figura 2. Sistema P-p para la columna se seccion variable 
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La matriz de rigidez se calculara mas adelante con Diferencias Finitas. 



1.2. Fundamento de las Diferencias Finitas 



Las diferencias finitas se fundamentan en encontrar expresiones 
aproximadas para las derivadas sean parciales o no. Se obtienen mediante el 
desarrollo de Taylor. Para el caso piano se tiene lo siguiente. Yepez (1985) 




Figura 3 Caso Piano de una Viga 
Se toman puntos simetricos todos equidistantes, como indica la figura 2. 
Si [/(q es una funcion diferenciable, al realizar el desarrollo de Taylor, se obtiene: 



Ui+i = Ui + h— + 



du h 2 d 2 u h 3 d 3 u 



+ 



d?yc 2 ! d^yc 3 ! d^yc 



+ 



U t -i= Ui-h— + 



du h 2 d 2 u h 3 d 3 u 



d^yc 2 ! d^yc 3 ! d^c 



+ 



Estas son las expresiones generales de las que se parte para hallar las 
aproximaciones de las derivadas en Diferencias Finitas. Para el exito de este 
metodo el paso h debe ser muy pequeno. Yepez (1985) 

Estas ecuaciones que aproximan las derivadas, se pueden representar en 
forma progresiva, regresiva o simetrica de acuerdo a la ubicacion del punto i 
(Punto Pivote). 

De esta forma el metodo consiste en la aproximacion diferencial en los 
n+2 puntos base (puntos que se toman dentro del problema) v 0 ,v 1 ,v 2 ,... v n ,v n+1 \ 
cada derivada tiene una representacion en diferencias finitas. De acuerdo a las 
condiciones de borde, se encuentra algunos valores de v t . Se tienen tantas 
condiciones de borde como grado de la ecuacion diferencial se quiera resolver. 
Una vez reemplazadas estas condiciones de borde en las ecuaciones en 
diferencias finitas, se conduce a un sistema de n ecuaciones simultaneas con n 
incognitas; la solucion de este sistema nos dara los valores restantes de las v, 
habiendo resuelto el problema. Yepez (1985) 

Mas adelante en el calculo de la matriz de rigidez de la viga se aclararan estos 
conceptos. 
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La nomenclatura a utilizar sera: 

du d 2 u d 3 u d 4 u 

—— = Ux , -— = Uxx , —j = Uxxx , -— - = Uxxxx 

LAsJL LAsJL LAssL LAssL 

Adoptando esta nomenclatura, se tiene: 

ft 2 ft 3 

U i+1 = Ui + h Ux + — Uxx + — Uxxx + ••• (1) 

/i 2 /i 3 

t/i-i = f/f — ft i/x + — f/xx — — Lfxxx + ••• (2) 

2. OBTENCION DE LA FORMULA PARA LA 1* DERIVADA EN 
DIFERENCIAS FINITAS. 

Como ejemplo para la investigacion, obtenemos las formulas progresiva, 
regresiva y simetrica, para la 1^ Derivada en Diferencias Finitas. 

2.1. Formula Progresiva 

Se trunca la ecuacion (1) en la primera derivada: 

U i+1 = Ui + hUx 

Ut+i ~ 

Ux = F. Progresiva 

ft 



2.2. Formula Regresiva 

Se trunca la ecuacion (2) en la primera derivada: 



= Ui-hUx 



U t - Ut-t 

Ux = F. Regresiva 

ft 



2.3. Formula Simetrica 

Se truncan las ecuaciones (1) y (2) en la segunda derivada: 

ft 2 

U i+1 = Ui + hUx + — Uxx (3) 
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h 2 

Ui-i = Ui — h Ux + — Uxx 



(4) 



Se resta (3) - (4): 



Ux = 



U i+1 - = 2hUx 
2h 



F. Simetrica 



2.4. Tablas de Derivadas en Diferencias Finitas 

Se presentan las tablas, hasta la 4^ Derivada en Diferencias Finitas, con 
sus respectivos errores. 

Tabla 1 Formulas Progresivas Ordinarias. Yepez (1985) 



Derivada 


Factor 


Ui 


u i+1 


U i+2 


U i+3 




Error 


Ux 


1 
h 


-1 


1 








— h Uxx 
2 


Uxx 


1 

h 2 


1 


-2 


1 






— h Uxxx 


Uxxx 


1 


-1 


3 


-3 


1 




3 /x \J ococococ 




h 3 








2 


Uxxxx 


1 

/i 4 


1 


-4 


6 


-4 


1 


—2h Uxx xxx 



Tabla 2 Formulas Regresivas Ordinarias. Yepez (1985) 



Derivada 


Factor 




Ui-3 


U t -2 


Ut-i 




Error 


Ux 


1 
h 








-1 


1 


h Uxx 
2 


Uxx 


1 

h 2 






1 


-2 


1 


h Uxxx 


Uxxx 


1 




-1 


3 


-3 


1 


3 h Uxxxx 




h 3 








2 


Uxxxx 


1 

/i 4 


1 


-4 


6 


-4 


1 


2 h Uxxxxx 
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Tabla 3 Formulas Simetricas Ordinarias. Yepez (1985) 



Derivada 


Factor 


U t -2 


Ut-t 




u i+1 


Ut-2 


Error 


Ux 


1 

2h 




-1 




1 




— h 2 Uxxx 
6 


Uxx 


1 

h 2 




1 


-2 


1 




— h 2 Uxxxx 
12 


Uxxx 


1 


-1 


2 




-2 


1 


— h 2 Uxxxxx 




2h 3 








4 


Uxxxx 


1 


1 


-4 


6 


-4 


1 


— h 2 Uxxxxx 




h 4 








6 


Se presentan tambien las tablas de Formulas Progresivas y Regresivas 
Mejoradas, que se utilizaran posteriormente en los calculos 




Tabla 4 Formulas Progresivas Me 


oradas. Yepez (1985) 


Derivada 


Factor 


u t 




U i+2 


U i+3 




Error 


Ux 


1 


-25 


48 


-36 


16 


-3 


I L K-J ✓v^v i/v i/v i/v 




12/i 




5 


Uxx 


1 


2 


-5 


4 


-1 




11 h 2 Uxxxx 




h 2 










12 


Uxxx 


1 


35 


-104 


114 


-56 


11 


—Sh 3 Uxxxxx 




12/i 2 


6 


7 I yyyy 


1 


-5 


18 


-24 


14 


-3 


7 /i 2 Ltoxxx 




2/i 3 






4 



Tabla 5 Formulas Regresivas Mejoradas. Yepez (1985) 



Derivada 


Factor 






Ut-2 


Ut-t 


u t 


Error 


Ux 


1 


-3 


16 


-36 


48 


-25 


4* r r /y- /y- /y- /y* /y* 


12/i 


5 


Uxx 


1 

/i 2 




1 


-4 


5 


-2 


11 h 2 Uxxxx 
12 
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Uxxx 


1 


11 


-56 


114 


-104 


35 


Sh 3 Uxxxxx 




12h 2 


6 


7 I / y 
K.J *\> *\> *\> *\> 


1 

2/i 3 


3 


-14 


24 


-18 


5 


7 h 2 Uxxxxx 
4 



2.5. Formulas Simetricas con Operadores 

Para esta investigacion se trabajara con las formulas simetricas ordinarias 
de las derivadas en Diferencias Finitas, para calculos mas exactos. 




2 h 



Uxx = 

h 2 




Uxxx = 



2 h 3 





T J 'yyyy ' — " 

LJ *\> *\> i/v i/v A 

3. ECUACION DIFERENCIAL PARA VIGAS DE SECCION 
VARIABLE SIN TOMAR EN CUENTA EL EFECTO DE CORTE 

La ecuacion diferencial para vigas de seccion variable es la siguiente. Yepez 
(1985): 

P ° = ~dx* + 2E ~dx~ ~dx^ + E ~dx^ ~dx~ 2 (5) 
Donde: 

I(x) ~* |n ^ r cia en el punto x (variable) 
v -> Desplazamiento vertical 
E -> Modulo de elasticidad del material 



P 0 -> Carga 



3.1. Planteamiento de la Ecuacion Diferencial de 4^ Orden con 
Diferencias Finitas 
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Para el calculo de la matriz de diferencias, se utiliza la 4^ derivada en Diferencias 
Finitas. 

Se toma un punto i (pivote) y tomando los operadores simetricos, antes 
vistos, se reemplaza en la ecuacion (5). 

n „. \ v i-2 ~ 4v f _i + 6v t - 4v i+1 + v i+2 i f-/i_i + h+i] \~Vi-2 + 2Vi-i - 2v i+1 + v i+2 ] 

Pi = EIi I £i \ + 2E [ — Th — Jl 2h* 

\ h-i ~ 2h + W l \ Vj-i ~ 2Vj + Vf+i ] 
+ I h 2 \[ h 2 \ 



Agrupando terminos en v, se tiene: 



Pi=E 



h-i + 2/f - /f 



i+l 



2/i 4 



-6/, + 2I i+1 



/i 4 
+ 2l t + 



Vi-i + E 



-2/i_! + 10/f - 2/ i+1 



/i 4 



Vf + E 



2/t-i - 6/ f - 



/i 4 



17; 



i+l 



2/i 4 



Sacando denominador comun 2/i 4 , se obtiene: 



2 P t K 



= Vt.Ah-i + 2/i - h +1 ] + v £ _ 1 [-12/( + 4/ £+1 ] + Vit-4/i.! + 20/* - 4/ £+1 ] + Pi +1 [4/;_ 1 - 12/J 
+ Vi_ 2 [-/i_ 1 + 2/i + / £+1 ] 



Con esta formula se obtienen los coeficientes de las incognitas del sistema 
de ecuaciones que se debe plantear para el calculo de la matriz de rigidez, para 
cualquier tipo de viga de seccion variable. 

Si los coeficientes son: 

B = -4/ £ _! + 20/, - 4I i+1 

C = M { _ x - 111, 

D = + 2/, + 

C = -12/ £ + 4/ £+1 

D' = / £ _! + 2/, - / £ + 1 

La ecuacion diferencial de 4^orden, quedaria de la siguiente forma. Yepez (1985): 



2 P t li 



= D'v t _ 2 + C' v ^ +Bv t + C v i+1 + D v t _ 2 



(6) 



4. CALCULO DE LA MATRIZ DE RIGIDEZ 



Para el calculo de la matriz de rigidez se recuerda que no pide carga, es decir: 
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P t = 0. 

Por lo tanto la ecuacion diferencial de 4^ orden en operadores, se presenta asi: 



E 




2h 4 

El paso h se debe tomar lo mas cercano a cero posible, para garantizar 
que la aproximacion del calculo sea lo mas exacta posible y es igual a la longitud 
de la viga analizada dividido para el numero de divisiones N. En la figura 3 se ha 
notado con i a un punto cualquiera en la viga, el que esta a la derecha sera el 
punto i + 1 y el que esta a la izquierda sera el punto i - 1. Yepez (1985) 



i-1 i i+1 




Figura 4 Discretizacion del dominio para la viga. 

El primer punto discreto de un elemento se tiene en el nudo inicial y a este 
se lo identifica con el numero 0 y el ultimo punto discreto se lo tiene en el nudo 
final identificado este como n. la ecuacion (6) se aplica desde el punto 1 hasta el 
punto n-1. No se aplica en el punto 0 ni en el punto n ya que son conocidos los 
giros en estos puntos. Aplicando la ecuacion (6) en el punto 1 se ve que se 
requiere de un punto auxiliar (punto ficticio) a la izquierda del punto 0, este punto 
no existe pero se necesita para poder aplicar la ecuacion (6). En la figura 4 se 
presenta el nudo inicial de la viga analizada identificado por el punto 0, el punto 
auxiliar a la izquierda identificado por -1, los puntos interiores con los que se va a 
trabajar, el punto en el nudo final identificado por 8 y el punto auxiliar a la derecha 
identificado por 9. Yepez (1985) 



En esta presentacion se han tornado 7 puntos discretos como se indica en la 
figura 4. 



-1 0123456789 




Figura 5 Identificacion de varios puntos discretos cercanos al nudo inicial. 

Mientras menos puntos auxiliares se consideren en la solucion del 
problema, se tendra mayor exactitud. Estos puntos discretos deben expresarse en 
funcion de los puntos reales utilizando para el efecto las condiciones de borde. 
Aguiar (2011) 
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La matriz de rigidez para el presente analisis quedana de la siguiente forma: 

_/c 2 i k 2 2. 

4.1. Calculo de la Primera Columna de k. 

Para el calculo de la primera columna de k, se deforma q 1 y se miden los 
momentos en los nudos inicial y final, respectivamente. 




Figura 6 Condiciones de borde para la primera columna de k. 



En la figura 5 se muestra las condiciones de borde tanto en el nudo inicial 
como en el nudo final, para la deformacion q 1 = 1. 

Se toma en cuenta que al partir de la ecuacion — = -— , se tiene la 
siguiente convencion de signos: 

• Giros horarios positivos. 

• Corrimiento vertical hacia abajo positive 

• Momento en el nudo inicial horario y momento en el nudo final anti 
horario positive 

• Corte en el nudo inicial hacia arriba y corte en el nudo final hacia abajo 
positive 

En el nudo inicial el giro 6 X=0 = -1, negativo por la convencion de signos. 
Si se toma la ecuacion del giro aproximada en operadores, se tiene: 




v -i = v ± + 2h 

En el nudo final por definicion el giro 6 X=L = 0, por lo que se obtiene: 
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Vn = Vc 



Una vez calculados los puntos auxiliares, se procede al calculo de la matriz de 
diferencias. 

• Cuando el punto pivote i = v 1 

D[ v_ 1 + B 1 v 1 + C 1 v 2 + D 1 v 3 = 0 
D[ (y ± + 2h) +B 1 v 1 + C ± v 2 + D 1 v 3 = 0 
{D[ + B x ) v 1 + C 1 v 2 + D 1 v 3 = -2 h D[ (7) 

• Cuando el punto pivote i = v 2 

C 2 v 1 +B 2 v 2 + C 2 v 3 +D 2 v 4 = 0 (8) 

• Cuando el punto pivote i = v 3 

D' 3 v t + C^v 2 + B 3 v 3 + C 3 v 4 + D 3 v 5 = 0 (9) 

• Cuando el punto pivote i = v 4 

D 4 v 2 + C 4 v 3 + B 4 v 4 + C 4 v 5 + D 4 v 6 = 0 (10) 

• Cuando el punto pivote i = v s 

Ds v 3 + v 4 + B 5 v 5 + C 5 v 6 +D 5 v 7 = 0 (11) 

• Cuando el punto pivote i = v 6 

D' 6 v 4 + C' 6 v s + B 6 v 6 + C 6 v 7 = 0 (12) 

• Cuando el punto pivote i = v 7 

D' 7 v s + C 7 v 6 + B 7 v 7 + D 7 v 9 = 0 
D' 7 v s + C 7 v 6 + (B 7 + D 7 ) v 7 = 0 (13) 
Finalmente se expresa el sistema de ecuaciones: 

S*v t = Q* (14) 
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c t 




0 


0 


0 


c 2 


B 2 


c 2 


D 2 


0 


0 


D' 3 


C 3 


B 3 


c 3 


D 3 


0 


0 


D\ 


C\ 


B 4 


Q 


D, 


0 


0 


D's 


C's 


B 5 


C 5 


0 


0 


0 


D'e 


C'e 


Be 


0 


0 


0 


0 


D> 7 


C 7 



0 
0 
0 
0 

c 6 

B 7 + D 7 \ 



-IV 




-2hD\- 


v 2 




0 


v 3 




0 


v 4 




0 


v s 




0 


v 6 




0 


-V 7 . 




0 



Donde S es la matriz de diferencias, v es el vector que contiene los 
desplazamientos verticales en cada punto discreto y Q* es el vector de cargas 
discretas. Se puede observar que la ecuacion (14) es similar a la ecuacion basica 
de estructuras, para analisis estatico K q = Q. Donde K es la matriz de rigidez, q 
es el vector de coordenadas y Q es el vector de cargas. 

Una vez resuelto el sistema se obtienen los v (desplazamientos 
verticales), con estos valores podemos calcular los valores de la primera columna 
de la matriz de rigidez. 



Asi se tiene que: 



k lt -> Es el momento en el nudo inicial. 
k 21 -> Es el momento en el nudo final 



d 2 v 

M = -El— 2 (15) 



Partiendo de la ecuacion (15) y respetando su convencion de signos, 
obtenemos valores de la primera columna de k 



fcn = E I 



o 



d 2 v 



(16) 



J 0 



/c 2 i — E I n 



d 2 v 
dx'' 



(20) 



Para obtener los valores de k xl y k 21 , con mayor precision, se trabaja con 
las ecuaciones de Derivadas en Diferencias Finitas mejoradas, progresivas para 
fen y regresivas para k 21 . 



4.2. Calculo de la Segunda Columna de k. 



Para el calculo de la segunda columna de k, se deforma q 2 y se miden los 
momentos en los nudos inicial y final, respectivamente. 
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Figura 7 Condiciones de borde para la segunda columna de k. 



En la figura 6 se muestra las condiciones de borde tanto en el nudo inicial 
como en el nudo final, para la deformacion q 2 = 1. 

En el nudo inicial por definicion el giro 6 X=0 = 0, por lo que se obtiene: 

v -i = v i 

En el nudo final el giro 6 X=L = -1, negativo por la convencion de signos. Si 
se toma la ecuacion del giro aproximada en operadores, se tiene: 




v 9 = v 7 — 2h 

Una vez calculados los puntos auxiliares, se procede al calculo de la matriz de 
diferencias. 

• Cuando el punto pivote i = v 1 

D[ v_ 1 + B 1 v 1 + C 1 v 2 + D 1 v 3 = 0 
(D[ + £ x ) v x + C 1 v 2 +D 1 v 3 = 0 (21) 

• Cuando el punto pivote i = v 2 

C 2 v 1 +B 2 v 2 + C 2 v 3 +D 2 v 4 = 0 (22) 

• Cuando el punto pivote i = v 3 

D 3 v 1 + C 3 v 2 + B 3 v 3 + C 3 v 4 + D 3 v 5 = 0 (23) 

• Cuando el punto pivote i = v 4 

D 4 v 2 + C 4 v 3 + B 4 v 4 + C 4 v 5 +D 4 v 6 = 0 (24) 
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Cuando el punto pivote i = v s 



Ds v 3 + C5 v 4 + £5 vs + C 5 v 6 +D 5 v 7 = 0 

Cuando el punto pivote i = v 6 
D' e v 4 + C' 6 v s + B 6 v 6 + C 6 v 7 = 0 

Cuando el punto pivote i = v 7 

D' 7 v s + C 7 v 6 + B 7 v 7 + D 7 v 9 = 0 

Z>7 v 5 + C7 i7 6 + (5 7 +D 7 )v 7 =2h D 7 



(25) 



(26) 



(27) 



Se expresa el sistema de ecuaciones: 



1 + B 1 




£>i 


0 


0 


c 2 




c 2 


D 2 


0 




c 3 


b 3 


c 3 




0 




c\ 


B 4 


Q 


0 


0 


D's 


C's 


s 5 


0 


0 


0 


D'e 


C'e 


0 


0 


0 


0 


D' 7 



0 
0 
0 

D 4 

Cs 
Be 
C 7 



0 
0 
0 
0 

c 6 









- 0 - 




v 2 




0 




V3 




0 




v 4 




0 




v s 




0 




v 6 




0 




lv 7 \ 




_2 h D 7 . 



Una vez resuelto el sistema se obtienen los v (desplazamientos 
verticales), con estos valores podemos calcular los valores de la segunda columna 
de la matriz de rigidez. 

Asi se tiene que: 

• k 12 -> Es el momento en el nudo inicial. 

• k 22 -> Es el momento en el nudo final 

Partiendo de la ecuacion (15) y respetando su convencion de signos, 
obtenemos valores de la segunda columna de k 



k 12 = -E I, 



o 



d 2 v 



(28) 



J 0 



k 22 = E I 



n 



d 2 v 



dx'' 



(29) 



Para obtener los valores de k 12 y k 22 , con mayor precision, se trabaja con 
las ecuaciones de Derivadas en Diferencias Finitas mejoradas, progresivas para 
^12 y regresivas para k 22 . 
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5. CODIGO EN MATLAB 

Se realizo un codigo en MATLAB, para resolver cada uno de los problemas 
expuestos en esta investigation. 

La viga de seccion variable analizada es la siguiente: 




Figura 8 Viga de seccion variable. 

Normalmente las cartelas en ambos extremos de la viga son de igual 
longitud (a = d). Para efecto de programacion del problema, se toma un caso lo 
mas general posible. 

El programa es amigable con el usuario, ya que se ingresan datos 
conocidos y calcula la matriz de rigidez directamente 

5.1. Ecuaciones de Altura e Inercia para la Viga de Seccion Variable 

Se presentan las ecuaciones en los tres tramos de la viga: 
• Tramo 1 : 0 < x < a 



m 0 = 



(hi - ho) 



a 



h x — h 0 
h t (x) = h 0 H * x 



a 



(30) 



/iW = 



12 



(31) 



Tramo 2: a < x < L - d 



h 2 (x) = /ii 



(32) 



h(x) = 



12 - 



(33) 
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Tramo 3: L - d < x < L 



K - K 



h 0 — h x 

h 3 (x) = — * x — h 

a 



o 



(34) 



b x * h 3 (x) 
12 



(35) 



5.2. Datos del Programa 



Los datos que se ingresan en el programa son los siguientes: 
• Las dimensiones de la seccion: 



o 


L 


— > 


Longitud de la viga. 


o 


a 




Longitud de la cartela a la izquierda. 


o 


c 


—> 


Longitud de la seccion constante. 


o 


d 


— > 


Longitud de la cartela a la derecha. 


o 


K 


— > 


Altura maxima en la seccion variable. 


o 


K 


—> 


Altura en la seccion constante. 


o 


b 


— > 


Espesor de la seccion. 



• El numero de puntos discretos en el dominio. 

o Este numero de puntos son solo los puntos que estan dentro 
del domino de la viga analizada, sin tomar en cuenta los 
puntos en los apoyos ni los puntos ficticios. 

• El modulo de elasticidad del material 



La columna de seccion variable analizada es la siguiente: 




Figura 9. Columna de seccion variable 
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5.3. Ecuaciones de Altura, Base e Inercia para la Columna de Seccion 
Variable 

• Variacion de la base 



b1 



(x) 



bo 




Figura 10. Variacion de la base en funcion de x (vista superior) 



bo — bl 



m = 



m m(x) m 

— = ; 772 (x) = — x 

L x L 



2m 2 (bo — bl) 
b(x) = bO — 2m(x); b(x) = bo — — x; b(x) = bo — — x 



(bo - bl) 
b{x) = bO x 



Variacion de la altura 
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L 


b1 










ho 


TaOO 






X 





Figura 11. Variacion de la altura en funcion de x (vista frontal) 



a a(x) a ho — hi 

- = ; a(x) = -X) a = 

L x L 2 

a 2 (ho — hi) 
h(x) = ho — 2a(x); /i(x) = ho — — x; /i(x) = ho — — x 

Li Li Li 

ho — hl\ 
1 x 



(no — ni\ 
/i(x) = ho — ( 1 



Calculo de la inercia 



/ = 



b * h 3 

b(x) * /i(x) 3 
12 



M -^^x)(/zo-(^i)x) 



3 
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5.4. Datos del programa 



E 
L 

hO 
hi 
bO 
b1 

pt 



Modulo de Elasticidad. 

Longitud Total de la viga. 

Altura para el punto inicial 

Altura para el punto final 

Ancho para el punto inicial del elemento 

Ancho para el punto final del elemento 

Numero de puntos discretos i 



Aislador sismico TADAS analizado es la siguiente: 




Figura 12. Modelo de un aislador sismico tipo TADAS 



5.5. Ecuaciones de Base e Inercia para el aislador sismico tipo TADAS 



• Calculo de la base 



bo — bi 



m = 



m 77i (x) m 

— = ; 771 (x) = —x 

L x L 



2m 2 (bo — bi) 
b(x) = bO — 2m(x); b(x) = bo — — x; b(x) = bo — x 



(bo — bi) 
b{x) = bO x 



• Calculo de la inercia 
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b*h 3 r ( b0 ~^ b ° l b ^ x ) h3 t b(x)*h 3 



1 12 1 12 1 12 

5.6. Datos del programa 

E Modulo de Elasticidad. 

h Altura del disipador 

bo Base mayor del disipador 

bi Base menor del disipador 

n Numero de placas 

t Espesor del disipador 

pt Numero de puntos discretos i 



6. EJERCICIO DE APLICACION 

Se desea calcular la matriz de rigidez de la siguiente viga, expuesta en la figura 




, 8.00 m 

Figura 13 Ejercicio de aplicacion. 

6.1. Metodo Anahtico 
6.1.1.Tablas de Guldan 



En el metodo de Guldan la matriz de rigidez se construye con las 
constantes de barra a, y b, para el problema expuesto se utilizara la Tabla 7 para 
cartelas rectas simetricas. 
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TabJa 7 

Carttlaj rccui limetricu. 

,4 



Constanta de barra a b 

X ==, Valor dc a (eacima) 



r 



BJ, 



Valor dc b (dcbajo) 0 « 



Lb 



Figura 13 Constantes de barra a y b. Guldan (1956) 

Donde / a es la inercia en la seccion constante de la viga, J A es la inercia 
en los extremos de la viga, l v es la longitud de las cartelas y / es la longitud total 
de la viga. 

60 * 80 3 
; = — _ — = 2560000 cm 4 

60 * 160 3 

J A = — = 20480000 cm 4 

JA 12 

J a 2560000 

n = — = = 0.125 

J A 20480000 

L 2.00 
X = JL= =0.25 
/ 8.00 

Con estos valores de n y A, ingresamos en la Tabla 6 (Tabla 7 de Guldan) 
y obtenemos las constantes de barra ay b. 

Tabla 6 Constantes de barra de Guldan para Cartelas Rectas Simetricas. Guldan 

(1956) 



A\n 


1,00 


0,90 


0,80 


0,70 


0,60 


0,50 


0,40 


0,30 


0,20 


0,15 


0,12 


0,50 


4,0 
2,0 


4,29 
2,19 


4,67 
2,43 


5,12 
2,73 


5,69 
3,12 


6,47 
3,66 


7,56 
4,43 


9,25 
5,68 


12,86 
8,04 


15,20 
10,25 


17,86 
12,48 


0,45 


4,0 
2,0 


4,28 
2,18 


4,63 
2,42 


5,05 
2,71 


5,59 
3,09 


6,31 
3,61 


7,81 
4,35 


8,84 
5,52 


11,58 
7,71 


14,04 
9,75 


16,31 
11,67 


0,40 


4,0 
2,0 


4,25 
2,18 


4,59 
2,41 


4,98 
2,68 


5,48 
3,04 


6,18 
3,53 


7,02 
4,21 


8,37 
5,28 


10,72 
7,22 


12,77 
8,97 


14,62 
10,57 


0,35 


4,0 
2,0 


4,25 
2,17 


4,54 
2,39 


4,90 
2,64 


5,35 
2,98 


5,92 
3,42 


6,70 
4,03 


7,86 
4,95 


9,79 
6,59 


11,41 

7,89 


12,83 
9,24 


0,30 


4,0 
2,0 


4,22 
2,16 


4,49 
2,35 


4,81 
2,59 


5,20 
2,89 


5,69 
3,27 


6,35 
3,80 


7,30 
4,58 


8,82 
5,87 


10,08 
6,92 


11,05 
7,83 


0,25 


4,0 
2,0 


4,20 
2,14 


4,43 
2,32 


4,70 
2,52 


5,03 
2,78 


5,44 
3,10 


5,97 
3,53 


6,71 
4,14 


7,84 
5,10 


8,69 
5,85 


9,38 
6,45 


0,20 


4,0 
2,0 


4,17 
2,13 


4,35 
2,27 


4,58 
2,44 


4,85 
2,65 


5,17 
2,91 


5,58 
3,24 


6,12 
3,89 


6,90 
4,35 


7,46 
4,84 


7,89 
5,21 


0,15 


4,0 


4,13 


4,28 


4,45 


4,65 


4,88 


5,17 


5,53 


6,08 


6,87 


6,62 
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2,0 


2,10 


2,22 


2,35 


2,51 


2,89 


2,93 


3,23 


3,65 


3,94 


4,16 


0,10 


4,0 
2,0 


4,09 
2,07 


4,20 
2,15 


4,51 
2,24 


4,44 
2,35 


4,59 
2,47 


4,76 
2,61 


4,98 
2,79 


5,26 
3,02 


5,44 
3,17 


5,57 
3,28 


0,05 


4,0 
2,0 


4,05 
2,04 


4,10 
2,08 


4,16 
2,13 


4,22 
2,18 


4,29 
2,23 


4,87 
2,30 


4,47 
2,37 


4,58 
2,47 


4,65 
2,53 


4,70 
2,57 


0,00 


4,0 
2,0 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 


4,00 
2,00 



Como vemos hay valores para X = 0.25, pero no existen para n = 0.125, 
Por lo que interpolamos los valores de n = 0.12 y n = 0.15. 

Tabla 7 Interpolacion de datos para ny A 



A/n 


0.12 


0.125 


0.15 




9.38 


9.265 


8.69 


0.25 










6.45 


6.35 


5.85 



Los valores obtenidos, haciendo interpolacion, para las constantes de 
barra son los siguientes: 

a = 9.265 
b = 6.350 

Una vez calculados los valores de las constantes de barras por el metodo 
de Guldan, se procede a armar la matriz de rigidez de la siguiente forma: 



K 



\a° b°] 



Donde: 



a 0 = 9.265 * 



b° = 6.35 * 



E*Ja 

L 

E*J a 



Los signos de cada valor de la matriz se consideran segun la misma ley de 
signos utilizada en todo el articulo. 



Resultados Metodo de las Tablas de Guldan 



Tomando como valor del modulo de elasticidad 1800000 T/m 2 , los 
resultados por el metodo analitico de Guldan son los siguientes: 
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It = 








1 


. 0e+Q4 * 






5 


.336640000000000 


-3 


.657600000000000 


-3 


. 657600000000000 


5 


.336640000000000 



Figura 14 Matriz de Rigidez por el Metodo de Guldan 
6.1.2.Tablas de Tena Colunga 

En el metodo de Tena Colunga es muy similar al metodo de Guldan donde 
la matriz de /rigidez se construye con los factores r llx , r 12x , r 22x ; que se toman de 
las tablas expuestas por el autor en el Apendice B del libro del que es autor. (Tena 
Colunga 2007). 

La tabla utilizada en este calculo fue la Tabla T26 del libro de Tena 
Colunga, de la cual obtenemos las constantes. 

Tabla 8 Resumen de Comparacion Tiempo y Error 




Tabla T26. Constantes de rigidez de elementos T de seccion variable, normalizados para Utk = 10, para a = 0.25 





Y 


r,u 




ru, 




r*. 




0.u 








025 


0.00 


3.79 


1.79 


3.79 


5.57 


5.57 


11.15 


24.00 


24.00 


12.00 


12.00 


0.20 


4.55 


2.33 


4.55 


6.8* 


6.88 


13.77 


25.11 


25.11 


11.34 


11.34 


0.40 


5.28 


2.89 


5.28 


8.17 


8.17 


16.34 


26.01 


26.01 


10.88 


10.88 


0.50 


5.63 


3.16 


5.63 


K.80 


8.80 


17.59 


26.41 


26.41 


10.69 


10.69 


0.60 


5.98 


3.44 


5.98 


9.41 


9.41 


18.82 


26.76 


26.76 


10.53 


10.53 


0.80 


6.63 


3.97 


6.63 


10.60 


10.60 


21 20 


27.40 


27.40 


10.28 


1028 


1.00 


7.25 


4.48 


7.25 


11.73 


11.73 


23 46 


27.94 


27.94 


10.07 


10.07 


1.25 


7.97 


5.08 


7.97 


13.05 


13.05 


26.10 


28.51 


28.51 


9.88 


9.88 


1 50 


863 


5.65 


8.63 


14.28 


14 28 


28.55 


29.00 


29.00 


9.73 


9.73 


1.75 


923 


6.17 


9.23 


15.41 


13.41 


30.81 


29.42 


29.42 


9.61 


9.61 


2.00 


9.79 


6.66 


9.79 


16.45 


16 45 


32.90 


29.78 


29.78 


9.51 


9.51 



De la Tabla 8 (Tabla T26 Tena Colunga), y con los datos del problema 
obtenemos los siguientes parametros: 

a = 0.25 



P = 0.25 
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y = 1.0 

L/ho = 10 

0.60 * 0.80 3 

L = — = 0.0256 m 4 

0 12 

Con estos daos preliminares, se procede a tomar los coeficientes de la 
Tabla 9 (Tabla T26 Tena Colunga). 

r iix ~ 7.25 
r 12x = 4.48 

r 22x ~ 

Con los parametros obtenidos, se arma la matriz de rigidez de la viga 
como esta expuesto en (Tena Colunga 2007). 




Resultados Metodo de las Tablas de Tena Colunga 

Tomando como valor del modulo de elasticidad 1800000 T/m 2 , los 
resultados por el metodo analftico de Tena Colunga son los siguientes, los signos 
son establecidos por la convencion de signos antes mencionada. 



k = 










1 


. Oe+04 * 








4 


.176000000000000 


-2 . 


530480000000000 




— Z 


.530480000000000 


4. 


17€000000000000 


1 



Figura 15 Matriz de Rigidez por el Metodo de Tena Colunga 
6.1. 3. Metodo de las Dovelas 

El metodo de las dovelas consiste en dividir a la viga analizada, en 
segmentos de areas rectangulares equivalentes, como lo indica la Figura 13. 
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Figura 16 Segmentacion de la Viga en Dovelas de Area Equivalente 

Se calculan las areas equivalentes de cada dovela rectangular y se da 
grados de libertad en cada nudo de las dovelas, como indica la Figura 13. 




Figura 17 Grados de Libertad de la Viga con Dovelas 

Se utiliza la condensacion de una matriz, expuesta en el capitulo 16 de 
Aguiar, 2014. La condensacion de una matriz consiste en colocar grados de 
libertad primarios o principales en la viga y grados de libertad secundarios, como 
se observa en la Figura 13. 

Al enumerar primero los grados de libertad primarios, se tiene en la matriz 
de rigidez total de la viga juntos a los valores que se van a necesitar una vez 
condensada la matriz, que en este caso sena de 6 x 6. 

Se utilizan los programas kmiembro y krigidez_por_elemento, de la lista de 
programas de CEINCI-LAB, para calcular la matriz k de rigidez con los 6 grados de 
libertad. Una vez calculada la matriz de rigidez se procede a condensar la matriz, 
para obtener la matriz de 2 x 2, la cual es la respuesta del problema planteado. 

La matriz condensada se consigue dividiendo la matriz de 6 x 6 en 
submatrices kaa, kba, kab y kbb, y se por ultimo se aplica la formula: 

^condensada ~ kaa — kab * kbb 1 * kba 

Se realizo un programa en MATLAB para calcular la matriz de rigidez por 
el metodo de las dovelas, en el cual se utilizaron lo programas de CEINCI-LAB 
antes mencionados. En el programa se ingresan los datos de cada seccion de las 
dovelas divididas, el vector de colocacion en cada segmento de la viga y su 
longitud. 
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Programa en MATLAB para el Metodo de las Dovelas 

%Programa para calcular la matriz de rigidez de una viga acartelada 
%por el metodo de las dovelas 
clear all 
clc 

ngl=6; E=1800000; seno=0; coseno=1; 

%Contribucion area 1 

b=0.60; h=1.40; L=1.0; 

[K3A1]=kmiembro(b,h 5 L 5 E,seno,coseno); 

VCA1=[0 0 1 0 0 3]; 

KT1=zeros(ngl); 

[KA1]=krigidez_por_elemento(KT1,K3A1 5 VCA1); 
%Contribucion area 2 
b=0.60; h=1.0; L=1.0; 
[K3A2]=kmiembro(b,h 5 L 5 E 5 seno,coseno); 
VCA2=[0 0 300 4]; 

[KA2]=krigidez_por_elemento(KA1,K3A2,VCA2); 
%Contribucion area 3 
b=0.60; h=0.80; L=4.0; 
[K3A3]=kmiembro(b,h 5 L 5 E,seno 5 coseno); 
VCA3=[0 0 400 5]; 

[KA3]=krigidez_por_elemento(KA2,K3A3,VCA3); 
%Contribucion area 4 
b=0.60; h=1.0; L=1.0; 
[K3A4]=kmiembro(b,h 5 L 5 E 5 seno,coseno); 
VCA4=[0 0 500 6]; 

[KA4]=krigidez_por_elemento(KA3,K3A4 5 VCA4); 
%Contribucion area 5 
b=0.60; h=1.40; L=1.0; 
[K3A5]=kmiembro(b,h 5 L 5 E,seno,coseno); 
VCA5=[0 0 600 2]; 

[KA5]=krigidez_por_elemento(KA4,K3A5,VCA5); 

%K TOTAL 

[KT]=KA5 

%CONDENSACION DE LA MATRIZ K 

na=2;kaa=KT(1 :na,1 :na);kab=KT(1 :na,na+1 :ngl);kba=KT(na+1 :ngl,1 :na); 
kbb=KT(na+1 :ngl,na+1 :ngl); 

k=kaa-kab*inv(kbb)*kba 

Resultados del Programa por el Metodo de las Dovelas 

Los resultados que se obtuvo con el metodo de las dovelas son los 
siguientes: 
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k = 








1.0e+05 








3.1658 


D 


.0002 




0.0002 


3 


.1658 





Figura 18 Resultados Metodo de las Dovelas 

Como se puede observar los resultados obtenidos no se comparan con los 
resultados en los anteriores metodos, por lo que el metodo de dovelas por efectos 
practicos, se descarta de la comparacion con respecto al metodo de Diferencias 
Finitas, ya que no se ajusta a los resultados. 

6.2. Codigo del Programa en Matlab 

Para el calculo de la matriz de rigidez por Diferencias Finitas se va a 
considerar 50 puntos discretos en la viga analizada para que la matriz de rigidez 
sea mas exacta y el modulo de elasticidad del material es 1800000 T/m 2 . El 
listado del programa que el usuario debe elaborar es el siguiente. 

% Programa para calcular la Matriz de Rigidez de una viga 
% de seccion variable doblemente acartelada por Dif. Finitas 
% con un P-p de 2 grados de libertad: Momento en cada uno 
% de los apoyos. Apoyos: Izquiera Fijo y Derecha Movil 

% 

% Por: Ronny Moreano 

% UFA - ESPE 

% 24 de mayo de 201 5 

% 

% Matriz_Rigidez_viga_doble_acartelada 
% 

format long 
tic 

%Datos del Problema 
E=1800000;%Modulo de Elasticidad (T/m2) 
L=8.00; %Longitud Total de la viga 
a=2.00; %Longitud de la cartela izquierda 
c=4.00; %Longitud de la parte continua de la viga 
d=2.00; %Longitud de la cartela derecha 
h0=1.60; %Altura en los extremos de la viga 
h1=0.80; %Altura en el centro de la viga 
b=0.60; %Espesor de la viga 
pt=50; %Numero de puntos discretos i 
N=pt+1 ; %IMumero de divisiones 
n=N+1; 

m0=(h1-h0)/a; %Pendiente en el tramo 1 
m1=(h0-h1)/d; %Pendiente en el tramo 3 
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%Paso h 
hini=L/N; 

%Calculo de las alturas e inercias en cada punto discreto 
for i=1 :n 

x(i)=(i-1)*hini; 

j=x(i); 
if j<=a 

h(i)=h0+m0*j; 
else if a<j && j<=L-d 
h(i)=h1; 
else j>L-d 

h(i)=m1*j-h0; 
end 
end 
end 

for i=1:n 

l(i)=(b*(h(i)) A 3)/12; 
end 

%Coeficientes de los Operadores de Dif. Finitas 
for i=2:pt+1 

B(i)=-4.*l(i-1 )+20.*l(i)-4.*l(i+1 ); 

C(i)=4*l(i-1)-12*l(i); 

D(i)=-l(i-1)+2*l(i)+l(i+1); 

C1(i)=-12*l(i)+4*l(i+1); 

D1(i)=l(i-1)+2*l(i)-l(i+1); 
end 

%Matriz de Diferencias 
S=zeros(pt); 
%Elemento S(1,1) 
for i=2 
for j=1 
for k=1 

S(j,k)=B(i)+D1(i); 
end 
end 
end 

%Diagonal Principal 
for j=2:pt-1 
for k=2:pt-1 
if j==k 

S(j,k)=B(j + 1); 
end 
end 
end 

%Diagonal Superior 1 
for j=1 :pt 
for k=2:pt 
if j==k-1 

S(j,k)=C(j + 1); 
end 
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end 
end 

%Diagonal Superior 2 
for j=1 :pt-2 
for k=3:pt 
if j==k-2 

S(j,k)=D(j + 1); 
end 
end 
end 

% %Diagonal Inferior 1 
for j=2:pt 
for k=1 :pt-1 
if j==k+1 

S(j,k)=C1(j+1); 
end 
end 
end 

% Diagonal Inferior 2 
for j=3:pt 
for k=1 :pt-1 
if j==k+2 

S(j,k)=D1(j + 1); 
end 
end 
end 

%Elemento S(pt,pt) 
for i=pt+1 
for j=pt 
for k=pt 

S(j,k)=B(i)+D(i); 
end 
end 
end 

S%Presenta la matriz de diferencias 
%Primera Columna de k 
M=zeros(1,pt); 
for i=2 

M(1)=-2*hini*D1(i); 

end 

Q1=M' 

v1=S\Q1%Desplazamientos en la viga 
a=zeros(2); 

a(1 ,1 )=E*I(1 )*(1 /((hini) A 2))*((2*v1 (1 ))-(5*v1 (2))+(4*v1 (3))-(v1 (4)));% Elemento 
(1,1) de la matriz de rigidez 

a(2,1 )=-E*l(n)*(1 /((hini) A 2))*((-2*v1 (pt))+(5*v1 (pt-1 ))-(4*v1(pt-2))+(v1(pt-3)));% 
Elemento (2,1) de la matriz de rigidez 

%Segunda Columna de k 

N=zeros(1,pt); 

for i=pt+1 
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N(pt)=2*hini*D(i); 

end 

Q2=N' 

v2=S\Q2%Desplazamientos en la viga 

a(1,2)=-E*l(1)*(1/((hini) A 2))*(^ Elemento 

(1,2) de la matriz de rigidez 

a(2,2)=E*l(n)*(1/((hini)^ 

Elemento (2,2) de la matriz de rigidez 

%Matriz de Rigidez k 
k=a 

Tiempo=toc 

Como se puede observar los unicos datos que se ingresan los 
anteriormente mencionados, esto se coloca en la seccion Datos del codigo. Una 
vez ingresados los datos, el programa calcula la matriz de rigidez directamente. 

6.2.1. Resultados del Programa 

El programa realiza los calculos y expone la matriz de rigidez k. Los 
resultados son los expuestos en la figura 15. 



Jc = 






1.0e+04 * 






5.276158043066232 


-3 


.716257011149362 


-3. 71625701115462 2 


5 


.276158043081701 


Tierapo = 






0.056895698631812 







Figura 19 Matriz de Rigidez k. 



Se desea calcular la matriz de rigidez de la siguiente columna, expuesta en la 
figura 20. 
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L 



Figura 20. Columna de seccion variable 

Datos: 

Modulo de Elasticidad 

E=21 00000 kg/cm A 2 

Longitud Total de la viga. 

L=6 m 

h0=0.6 m 

h1=0.3 m 

b0=0.60 m 

b1=0.30 m 

Numero de puntos discretos 6 



6.3. Metodo de las diferencias finitas 

6.3.1. Programa realizado en MATLAB 

% UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE 

% Calculo la Matriz de Rigidez k para un elemento de seccion variable tanto 

% en la base como la altura 

% Autor: Dennisse Criollo. 

% Estudiante de V nivel de Ingenieria Civil. 

% Departamento de Ciencias de La tierra y la Construccion. 

% Julio 2015. 

format long 

%Datos del Problema. 

clear all 

clc 

E=21 00000; %Modulo de Elasticidad. 
L=5; %Longitud Total de la viga. 

h0=0.60; %Altura para el punto inicial . 
h1=0.30; %Altura para EL punto final. 
b0=0.60; %Ancho para el punto incial del elemento 
b1 =0.30; %Ancho para el punto final del elemento 
pt=5; %Numero de puntos discretos i 

N=pt+1 ; %Numero de divisiones 
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n=N; 
P=2; 

q1=(h0-h1)/2; 
hipo=sqrt(q1 A 2+L A 2); 

%Paso h 
hini=L/N; 

%Calculo de las alturas e inercias en cada punto discreto 

for i=1 :pt+2 

pass=(i-1)*hini; 
b(i)= bO-(pass*((bO-b1)/L)); 
H(i)= hO-(pass*((hO-h1)/L)); 
l(i)=(b(i)*(H(i)) A 3)/12; 
end 
I; 

%Coeficientes de los Operadores de Diferencias Finitas 
for i=2:pt+1 

B(i)=-4.*l(i-1)+20.*l(i)-4.*l(i+1); 

C(i)=4*l(i-1)-12*l(i); 

D(i)=-l(i-1)+2*l(i)+l(i+1); 

C1(i)=-12*l(i)+4*l(i+1); 

D1(i)=l(i-1)+2*l(i)-l(i+1); 
end 

%Matriz de Diferencias 
S=zeros(pt); 
%Elemento S(1,1) 
for i=2 
for j=1 
for k=1 

S(j 3 k)=B(i)+D1(i); 
end 
end 
end 

%Diagonal Principal 
for j=2:pt-1 
for k=2:pt-1 
if j==k 

S(j,k)=B(j + 1); 
end 
end 
end 

%Diagonal Superior 1 

for j=1 :pt-1 
for k=2:pt 
if j==k-1 

S(j,k)=C(j + 1); 
end 
end 
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end 

%Diagonal Superior 2 
for j=1:pt-2 
for k=3:pt 
if j==k-2 

S(j,k)=D(j + 1); 
end 
end 
end 

%Diagonal Inferior 1 
for j=2:pt 
for k=1 :pt-1 
if j==k+1 

S(j,k)=C1(j + 1); 
end 
end 
end 

% Diagonal Inferior 2 
for j=3:pt 
for k=1 :pt-1 
if j==k+2 

S(j,k)=D1(j + 1); 
end 
end 
end 

%Elemento S(pt,pt) 
for i=pt+1 
for j=pt 
for k=pt 

S(j 3 k)=B(i)+D(i); 
end 
end 
end 

S ;%Presenta la matriz de diferencias 
Q=(2*P*((hini)M))/E; 
P=zeros(pt,1); 
for j=1 :pt 

P(jM)=Q; 
end 

v=GAUSSJORD(S,P); 

% %Primera Columna de k 
M=zeros(1 ,pt); 
for i=2 

M(1)=-2*hini*D1(i); 

end 

Q1=M'; 

v1 =S\Q1 ;%Desplazamientos en la viga 
%Segunda Columna de k 
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N=zeros(1 ,pt); 

for i=pt+1 

N(pt)=2*hini*D(i); 

end 

Q2=N'; 

v2=S\Q2 ;%Desplazamientos en la viga 

%MATRIZ DE RIGIDEZ k 
k=zeros(2,2); 

a(1 ,1 )=(E*I(1 )*(-1 04*v1 (1 )+1 1 4*v1 (2)-56*v1 (3)+1 1 *v1 (4)))/(1 2*hini A 2); 

a(2,1 )=(E*I(1 )*(-1 04*v2(1 )+1 1 4*v2(2)-56*v2(3)+1 1 *v2(4)))/(1 2*hini A 2); 

% a(1 ,2)=(E*I(1 )*(-1 04*v2(1 )+1 1 4*v2(2)-56*v2(3)+1 1 *v2(4)))/(1 2*hini A 2); 

%a(2,2)=-(E*l(1 04*v2(pt)+1 1 4*v2(pt-1 )-56*v2(pt-2)+1 1 *v2(pt-3)))/(1 2*hini A 2) ; 

a(1,2)=a(2,1) 

a(2,2)=a(1,1) 

k=a 

%grafico de la viga 
subplot(3,1,2); 
%axis([0 L+2-1 hO+1]) 
x1=0:hini:L; 
x2=[0 L]; 
y2=[h0q1+h1]; 
y1=((hO-h1)/(2*L))*x1; 
hold on 
f=[LL]; 
g=[q1 h0-q1 ]; 

plot (x1 ,y1 5 x2 5 y2;b',f ,g/b7LineWidth\2) 
grid on 

ylabelCAItura 1 ) 

title('Cara frontal de la columna de seccion variable ') 
xlabel('Longitud') 

%GRAFICO DE DESPLAZAMIENTOS 
subplot(3,1,3); 
for t=1 :pt 

u(t)=hini*t; 
end 
hold on 
m=[0 hini]; 
n=[0 v(1)]; 
p=[hipo hipo-hini]; 
q=[0 v(pt)]; 

plot^rV/r'^rn/r'prq/r'/LineWidth'^) 
grid on 

ylabel('Desplazamientos') 
title('Graficas de deformacion ') 
legend('Desplazamientos con 15 puntos ',-1) 
xlabel('Longitud') 

%Comparacion con 100 pts 
v4=v8(100); 
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pt1=100; 
pass=L/(pt1+1); 
for r2=1 :pt1 

u2(r2)=pass*r2; 
end 
hold on 
t1=[0 pass]; 
r2=[0 -v4(1)]; 
f3=[L-pass L]; 
f4=[-v4(pt1) 0] ; 
subplot(3,1,3); 
plot(u2,-v4/b') 
%Comparacion con 200 pts 
v6=v8(200); 
pt2=200; 
pass=L/(pt2+1); 
for r2=1 :pt2 

u3(r2)=pass*r2; 
end 
hold on 
t5=[0 pass]; 
r6=[0 -v4(1)]; 
f7=[L-pass L]; 
f8=[-v4(pt1) 0] ; 
subplot(3,1,3); 
plot(u3,-v6,'g') 

%IMAGEN 

k1=imread(Viga1 .jpg'); 
subplot(3,1 ,1 ) 
imshow(k1,[0 500]) 



6.3.1 .1 . Resultados del programa 



1.0e+003 * 

7.644448306582492 1.185359111378156 
1.185359111378156 7.644448306582492 



Figura 21. Matriz de rigidez 
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Figura 22. Imagen, grafico de la vista frontal en 2D, y desplazamientos de la 

columna de seccion variable. 



Se desea calcular la matriz de rigidez del siguiente aislador sfsmico tipo TADAS 
expuesto en la figura 22 -23. 

y bo = 0.30 m l 



o 

CD 




pi = 0.15 my 

Figura 22. Modelo aislador sfsmico tipo TADAS posicion vertical 
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Figura 22. Modelo aislador sismico tipo TADAS posicion horizontal 

6.4. Solucion en diferencias finitas 

Se divide al elemento en puntos discretos, para el ejemplo de aplicacion 
se utilizo 6 puntos de discretizacion. 





Figura 23. Modelo aislador sismico tipo TADAS con puntos discretos 



Se realiza el procedimiento de diferencias finitas para hallar la matriz de 
rigidez con ayuda del programa realizado en MATLAB 

6.4.1. Codigo del programa en MATLAB 

% UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE 

% Calculo la Matriz de Rigidez k para un elemento de seccion variable- 

% aplicacion aislador sisimico tipo TADAS 

% Autor: Dennisse Criollo. 

% Estudiante de V nivel de Ingenieria Civil 

% Departamento de Ciencias de La tierra y la Construccion. 

% Julio 2015. 
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format long 

%Datos del Problema. 

clear all 

clc 

E=2100000; %Modulo de Elasticidad. 

L=0.60; %Longitud Total de la viga. L=h (m) 

b0=0.30; %Base en el punto inicial del elemento (m) 

bl=0.15; %Base en el punto final del elemento. bl=bi (m) 

t=0.03; %Espesor constante del elemento (m) 

pl=6; % Numero de placas 

pt=75; %Numero de puntos discretos i 

N=pt+1; %Numero de divisiones 

n=N; 

P=2; 

ql=(b0-bl)/2; 
hipo=sqrt(ql A 2+L A 2); 
%Paso h 

hini=L/N; %hini =k 
p=hini; 

%Calculo de las alturas e inercias en cada punto discreto 

for i=l:pt+2 

pass=(i-l)*hini; 
b(i)= bO-(((bO-bl)/L)*pass); 

H(i)=t*pl; % es igual al espesor t*pl (numero de placas) 

l(i)=(b(i)*(H(i)) A 3)/12; 

end 

I; 

%Coeficientes de los Operadores de Diferencias Finitas 
for i=2:pt+l 

B(i)=-4.*l(i-l)+20.*l(i)-4.*l(i+l); 

C(i)=4*l(i-l)-12*l(i); 

D(i)=-l(i-l)+2*l(i)+l(i+l); 

Cl(i)=-12*l(i)+4*l(i+l); 

Dl(i)=l(i-l)+2*l(i)-l(i+l); 
end 

%Matriz de Diferencias 
S=zeros(pt); 
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%Elemento S(l,l) 
for i=2 
forj=l 
for k=l 

S(j,k)=B(i)+Dl(i); 
end 
end 
end 

%Diagonal Principal 
for j=2:pt-l 
for k=2:pt-l 
if j==k 

S(j,k)=B(j+l); 
end 
end 
end 

%Diagonal Superior 1 

for j=l:pt-l 
for k=2:pt 
if j==k-l 

S(j,k)=C(j+l); 
end 
end 
end 

%Diagonal Superior 2 
for j=l:pt-2 
for k=3:pt 
if j==k-2 

S(j,k)=D(j+l); 
end 
end 
end 

%Diagonal Inferior 1 
for j=2:pt 
for k=l:pt-l 
if j==k+l 

S(j,k)=Cl(j+l); 
end 
end 
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end 

% Diagonal Inferior 2 
for j=3:pt 
for k=l:pt-l 
if j==k+2 

S0,k)=Dl(j+l); 
end 

end 

end 

%Elemento S(pt,pt) 
for i=pt+l 
for j=pt 
for k=pt 

S(j,k)=B(i)+D(i); 
end 
end 
end 

S; %Presenta la matriz de diferencias 
Q=(2*P*((hini) A 4))/E; 

P=zeros(pt,l); 
for j=l:pt 

P(j,l)=Q; 
end 

v=GAUSSJORD(S,P); 

%MATRIZ DE RIGIDEZ k 
%Vectores de cargas verticales 
Qi=zeros(pt,l); 

Qi(l / l)=-2*p*(l(l)+2*l(2)-l(3)); 
Qf=zeros(pt,l); 

Qf(pt # l)=2*p*(-l(n-2)+2*l(n-l)+l(n)); 

%Vectores de desplazamientos verticales 

vi=inv(S)*Qi; 

vf=inv(S)*Qf; 

%Matriz de rigidez 
k=zeros(2,2); 



Calculo de la matriz de rigidez de un elemento ... 303 



a(l,l)=E*l(l)*(l/((hini) A 2))*((2*vi(l))-(5*vi(2))+(4*vi(3))-(vi(4)));% 
Elemento (1,1) de la matriz de rigidez 

a(2,l)=E*l(n)*(l/((hini) A 2))*((-2*vi(pt))+(5*vi(pt-l))-(4*vi(pt-2))+(vi(pt- 
3)));% Elemento (2,1) de la matriz de rigidez 

a(l,2)=E*l(l)*(l/((hini) A 2))*((2*vf(l))-(5*vf(2))+(4*vf(3))-(vf(4)));% 
Elemento (1,2) de la matriz de rigidez 

a(2,2)=E*l(n)*(l/((hini) A 2))*((-2*vf(pt))+(5*vf(pt-l))-(4*vf(pt-2))+(vf(pt- 

3)));% Elemento (2,2) de la matriz de rigidez 

k= 

6.4.2. Resultados del programa 

k = 

1.0e+003 * 

1.729588983718263 0.718317786118304 
0.696280790636672 1.229259675713808 

Figura 24. Matriz de rigidez metodo de diferencias finitas 

6.5. Solucion metodo de dovelas para el modelo de aislador si'smico tipo 
TADAS 

v- 

nt =6*0.03 =0.18 




Figura 25. Esquema de la division de dovelas para el modelo aislador sfsmico 

tipo TADAS 
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Figura 26. Modelo aislador sismico tipo TADAS con grados de libertad principales 

y secundarios 

6.5.1. Codigo del programa en MATLAB para dovelas 

%UNIVERSIDAD DE LAS FUERZAS ARMADAS - ESPE 
%DEPARTAMENTO DE CIENCIAS DE LA TIERRA Y LA CONSTRUCCION 
%INGENIERIA CIVIL 

%ANALISIS MATRICIAL DE ESTRUCTURAS 

%METODO DE DOVELAS PARA ELEMENTO DE SECCION VARIABLE 

clc 

clear 

format short 

%DATOS 

Kt=zeros(18); 

B=0.3; %Ancho de la viga 

LT=0.6; %Longitud total de la viga 

L=LT/5; %Longitud de cada dovela cte 

nt=0.18; %altura = n*t;donde n es el numero de placas y t = espesor 
E=2100000; 

%MIEMBRO 1 
coseno=l; 
seno=0; 
bl=0.1*B; 

[K3]=kmiembro(bl,nt,L,E,seno,coseno); 
VC=[1 2 3 7 8 9]; 

[Ktl]=krigidez_por_elemento(Kt,K3,VC); 
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%MIEMBRO 2 
coseno=l; 
seno=0; 
b2=0.3*B; 

[K31]=kmiembro(b2,nt,L,E,seno,coseno); 
VC=[7 8 9 10 11 12]; 

[Kt2]=krigidez_por_elemento(Ktl,K31,VC); 

%MIEMBRO 3 
coseno=l; 
seno=0; 
b3=0.5*B; 

[K32]=kmiembro(b3,nt,L,E,seno,coseno); 
VC=[10 11 12 13 14 15]; 
[Kt3]=krigidez_por_elemento(Kt2,K32,VC); 

%MIEMBRO 4 
coseno=l; 
seno=0; 
b4=0.7*B; 

[K33]=kmiembro(b4,nt,L,E,seno,coseno); 
VC=[13 14 15 16 17 18 ]; 
[Kt4]=krigidez_por_elemento(Kt3,K33,VC); 

%MIEMBRO 5 
coseno=l; 
seno=0; 
b5=0.9*B; 

[K34]=kmiembro(b5,nt,L,E,seno,coseno); 
VC=[16 17 18 4 5 6]; 

[KT]=krigidez_por_elemento(Kt4,K34,VC) 

%Condensaci6n de la matriz de rigidez a 
kaa=KT(l:6,l:6); 

kab=KT(l:6,7:18); 
kba=kab'; 

kbb=KT(7:18,7:18); 
KL=kaa-kab*inv(kbb)*kba 



306 



Moreano Ronny , Criollo Sandoval Dennisse 



6.5.2. Resultados del programa 

KL = 

1.0e+004 * 

0 0 

0.4535 0.1960 

0.0761 0.0186 

0 0 

0.4535 -0.1960 

0.1960 0.0990 



5.2873 0 0 -5.2873 

0 0.4535 0.0761 0 

0 0.0761 0.0271 0 

-5.2873 0 0 5.2873 

0 -0.4535 -0.0761 0 

0 0.1960 0.0186 0 



Figura 27. Matriz de rigidez condensada lateral metodo de dovelas 

7. COMPARACION DE METODOS 

Se basara la comparacion de metodos con la viga acartelada de seccion 
variable con el fin de determinar cual de ellos es el mas factible. 

Se tiene como resultado analitico el valor exacto de la matriz de rigidez, 
para la viga acartelada y se sabe que los valores obtenidos con el metodo de 
diferencias finitas dependen del numero de puntos discretos que se asigne en el 
programa. 

A continuacion se muestra el analisis de comparacion, donde se revisa el 
error del metodo numerico (Diferencias Finitas) asignando distintos valores al 
numero de puntos discretos, con respecto a la respuesta del metodo analitico 
(Metodo de Guldan). 

7.1. Resultados del Programa de Matlab 

7.1 .1 . Matriz de Rigidez para 7 Puntos Discretos 



k: = 








1.0e+05 * 








6.343366948601166 


-4. 


809368900455437 




-4.80936890045543S 


6. 


348366948601181 




Tiempo = 








0.067296199155900 









Figura 28 Matriz de Rigidez para 7 Puntos Discretos 
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7.1.2. Matriz de Rigidez para 25 Puntos Discretos 

k = 

l.Qe+04 * 

5 . 297132434597978 -3 . 720331400943342 
-3 . 720331400943401 5 . 297132434598152 

Tiempo = 

0.059476579415648 

Figura 29 Matriz de Rigidez para 25 Puntos Discretos 

7.1.3. Matriz de Rigidez para 50 Puntos Discretos 



k = 






1.0e+04 * 






5.276158043066232 


-3 


.716257011149362 


-3.716257011154622 


5 


.276158043081701 


Tierapo = 






0.056895698631812 







Figura 30 Matriz de Rigidez para 50 Puntos Discretos 
7.1.4. Matriz de Rigidez para 100 Puntos Discretos 

k = 

1.0e+04 * 

5.313373556966717 
-3 . 695250496631249 

Tieirpo = 

0 . 121535750223769 



-3.695250496631051 
5 . 313873556968348 



Figura 31 Matriz de Rigidez para 100 Puntos Discretos 
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7.1.5. Matriz de Rigidez para 150 Puntos Discretos 



k = 






1.0e+04 * 






5.317858745450432 
-3. 634566898907217 


-3. 
5. 


684566898833754 
317858745011180 


Tierapo = 






0.337762669302144 







Figura 32 Matriz de Rigidez para 150 Puntos Discretos 
7.1.6. Matriz de Rigidez para 1000 Puntos Discretos 

k = 

1.0e+04 * 

5.325642349789341 
-3 . 656929734127979 

Tiercpo = 

8.511961217080941 

Figura 33 Matriz de Rigidez para 1000 Puntos Discretos 
Matriz de Rigidez para 3000 Puntos Discretos 

k = 

1.0e+04 * 

5 . 326311613459200 -3 . 653373161245704 

-3 . 653373649583176 5 . 32 62 3055 9 994335 

T i emp □ = 

31 . 358899804947350 

Figura 34 Matriz de Rigidez para 3000 Puntos Discretos 



-3. 65692 9738011689 
5.325642474110362 
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Matriz de Rigidez para 5000 Puntos Discretos 



k = 






1.0e+04 * 






5.325039955444503 
-3 .653748783189208 


-3. 
5. 


653747179308931 
325035038972200 


Tieir.po = 






78 .668016621984137 







Figura 35 Matriz de Rigidez para 5000 Puntos Discretos 



7.2. Resumen y Analisis de Comparacion 

Primero se comparan los dos metodo analiticos (Guldan vs Tena 
Colunga), donde podemos apreciar ya un error. Si comparamos los coeficientes de 
la matriz y calculamos el error absolute, tenemos que: 

\x -x A | 
E = - — — — - = 0.25599 
\x\ 

Esto ya representa un error pequeno donde se tendra que observar los 
resultados del metodo de las diferencias finita, para poder deducir a cual de los 
dos metodos se ajusta mas Diferencias Finitas. 

A continuacion se realiza el calculo del error con los metodos analiticos 
propuestos anteriormente. 

7.2.1. Analisis del Metodo de Diferencias Finitas vs Guldan 

Tabla 9 Resumen de Comparacion Tiempo y Error 



# Puntos 
Discretos 


Tiempo 


Error 


7 


0,06729 


0,21453 


25 


0,05948 


0,01022 


50 


0,05690 


0,01176 


100 


0,12154 


0,00605 


150 


0,33776 


0,00456 


1000 


8,51196 


0,00110 


3000 


31,35890 


0,00141 


5000 


78,66802 


0,00149 



Se puede observar como el error va siendo menor cuantos mas puntos 
discretos se le asignen al programa. 
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Con los resultados obtenidos, se generan graficos de los tiempos y errores 
de calculo para el metodo de las Diferencias Finitas, obteniendo asi las tendencias 
en las que se desenvuelve cada variable. 



E 



# Puntos Discretos - Tiempo de Calculo en 

Matlab 

90.00000 | 1 1 1 1 1 — 

80.00000 1— 

70.00000 ~S\ — 

60.00000 

50.00000 

40.00000 

30.00000 

20.00000 

10.00000 • 

0.00000 mm 
-10.00000 0 1000 



2000 3000 4000 
# Puntos Discretos 



5000 



6000 



Figura 36 Grafico: # Puntos Discretos - Tiempo 

# Puntos Discretos - Error de Calculo en 

Matlab 
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Figura 37 Grafico: # Puntos Discretos - Error 



En las figuras 24 y 25, se muestran la tendencia de cada variable del 
metodo de las Diferencias Finitas con respecto al metodo analitico; se puede 
verificar que mientras mas puntos discretos se asignen en el programa el error va 
disminuyendo, pero el tiempo de calculo aumenta, esto se debe a que el programa 
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de Matlab realiza mas operaciones y mientras mas puntos discretos se asignen 
mas se va a tardar en calcular las respuesta. 

Por esta razon se puede decir que si bien mientras mas puntos discretos 
se asignen en el programa se genera menos error, pero tambien el programa se 
tarda mas en dar su resultado; por lo que no se debe asignar ni pocos ni muchos 
puntos discretos. 

Observando los resultados del analisis, se busca un numero de puntos 
discretos que no genere demasiado error y que el tiempo de calculo no sea alto, 
para optimizar el proceso; este tiempo tambien puede depender del ordenador con 
el que el usuario trabaje. 

Como minimo se asigna 50 puntos discretos al proceso, pero se puede 
tomar como media 1000 puntos, para optimizar el calculo y que la respuesta se 
acerque mas a la exacta. 

7.2.2. Analisis de Metodo de Diferencias Finitas vs Tena 

Tabla 10 Resumen de Comparacion Tiempo y Error 



# Puntos 
Discretos 


Tiempo 


Error 


7 


0,06729 


0,26015 


25 


0,05948 


0,13425 


50 


0,05690 


0,13175 


100 


0,12154 


0,13626 


150 


0,33776 


0,13674 


1000 


8,51196 


0,13766 


3000 


31,35890 


0,13775 


5000 


78,66802 


0,13759 



Como se observa el tiempo de calculo del programa en Matlab no cambia, 
lo que cambia con relacion al anterior metodo es el analisis del error. Se observa 
en el error que la matriz de rigidez obtenida por Diferencias Finitas genera con 
respecto al metodo de Tena (2007), se mantiene constante en la mayoria de la 
discretizacion por puntos. Esto genera el siguiente grafico de analisis. 



312 



Moreano Ronny , Criollo Sandoval Dennisse 



# Puntos Discretos - Error de Calculo en 

Matlab 

0.30000 |— 
0.25000 # 
0.20000 I— 
§ 0.15000 ^ 
0.10000 U 
0.05000 U 
0.00000 — 

0 1000 2000 3000 4000 5000 6000 

# Puntos Discretos 

Figura 38 Grafico: # Puntos Discretos - Error 

Como se verifica en la Figura 26, podemos ver que el error se mantiene 
constante casi en un valor. 

Por esta razon podemos concluir que el metodo de Tena no se ajusta tan 
bien como el metodo de Guldan. 

8. COMENTARIOS Y CONCLUSIONES 

Se ha presentado la deduccion de alguna de las formulas y operadores en 
Diferencias Finitas, con el fin de dar al lector un acercamiento mas profundo a la 
teoria de este metodo. Se observa que el metodo no es demasiado complicado y 
puede ser dominado. 

En este artfculo se procede con el calculo de la matriz de rigidez para un 
tipo de viga de seccion variable, columna de seccion variable y elementos 
aisladores sismicos tipo TADAS que se presenta con frecuencia en la Ingeniena 
Civil. 

Se resuelve el problema analiticamente por los metodos de Guldan, Tena- 
Colunga, y dovelas, y se observa que el metodo analitico de las Tablas de Guldan 
es el mas exacto para el caso de la viga de seccion variable. Mientras que para la 
columna de seccion variable y el aislador sismico tipo TADAS se realizo la 
comparacion unicamente con el metodo de dovelas observando que es mucho 
mas preciso el metodo por diferencias finitas. 

Se ha presentado el codigo en Matlab, y se ha detallado la manera de 
ingresar los datos. Se puede ver que el programa es bastante sencillo de usar y de 
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ingresar los datos del problema. Con este programa se aspira aportar a la 
ensenanza de los metodos de Diferencias Finitas aplicados a la Ingeniera Civil. 

Finalmente se realiza una comparacion de los dos metodos (analitico y 
numerico), con el fin de comprobar los valores de la matriz de rigidez y la 
funcionalidad y eficacia del metodo numerico de las Diferencias Finitas. Se 
comprueba que el metodo de Diferencias Finitas, primero depende del numero de 
puntos discretos que se le asigne, y segundo al dar un numero de puntos discretos 
alto, el metodo responde favorablemente con un error casi imperceptible. 

Se recomienda como conclusion, en el metodo de Diferencias Finitas, 
trabajar con 1000 puntos discretos, para obtener un resultado favorable para el 
calculo de la matriz de rigidez, y con un tiempo de calculo del programa no tan 
grande y optimizado. 
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